﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Data;

namespace BaiduMapArea
{
    class GetArea
    {
        public static Double GetPolygonArea(DataTable pts)
        {
            double EARTHRADIUS = 6370996.81;
            double totalArea = 0;//初始化总面积
            double LowX = 0.0;
            double LowY = 0.0;
            double MiddleX = 0.0;
            double MiddleY = 0.0;
            double HighX = 0.0;
            double HighY = 0.0;
            double AM = 0.0;
            double BM = 0.0;
            double CM = 0.0;
            double AL = 0.0;
            double BL = 0.0;
            double CL = 0.0;
            double AH = 0.0;
            double BH = 0.0;
            double CH = 0.0;
            double CoefficientL = 0.0;
            double CoefficientH = 0.0;
            double ALtangent = 0.0;
            double BLtangent = 0.0;
            double CLtangent = 0.0;
            double AHtangent = 0.0;
            double BHtangent = 0.0;
            double CHtangent = 0.0;
            double ANormalLine = 0.0;
            double BNormalLine = 0.0;
            double CNormalLine = 0.0;
            double OrientationValue = 0.0;
            double AngleCos = 0.0;
            double Sum1 = 0.0;
            double Sum2 = 0.0;
            double Count2 = 0;
            double Count1 = 0;
            double Sum = 0.0;
            double Radius = EARTHRADIUS; //6378137.0f,WGS84椭球半径 
            int Count = pts.Rows.Count;
            
            for (int i = 0; i < Count; i++)
            {
                if (i == 0)
                {
                    LowX = Convert.ToDouble(pts.Rows[Count - 1][0].ToString()) * Math.PI / 180;
                    LowY = Convert.ToDouble(pts.Rows[Count - 1][1].ToString()) * Math.PI / 180;
                    MiddleX = Convert.ToDouble(pts.Rows[0][0].ToString()) * Math.PI / 180;
                    MiddleY = Convert.ToDouble(pts.Rows[0][1].ToString()) * Math.PI / 180;
                    HighX = Convert.ToDouble(pts.Rows[1][0].ToString()) * Math.PI / 180;
                    HighY = Convert.ToDouble(pts.Rows[1][1].ToString()) * Math.PI / 180;
                }
                else if (i == Count - 1)
                {
                    LowX = Convert.ToDouble(pts.Rows[Count - 2][0].ToString()) * Math.PI / 180;
                    LowY = Convert.ToDouble(pts.Rows[Count - 2][1].ToString()) * Math.PI / 180;
                    MiddleX = Convert.ToDouble(pts.Rows[Count - 1][0].ToString()) * Math.PI / 180;
                    MiddleY = Convert.ToDouble(pts.Rows[Count - 1][1].ToString()) * Math.PI / 180;
                    HighX = Convert.ToDouble(pts.Rows[0][0].ToString()) * Math.PI / 180;
                    HighY = Convert.ToDouble(pts.Rows[0][1].ToString()) * Math.PI / 180;
                }
                else
                {
                    LowX = Convert.ToDouble(pts.Rows[i - 1][0].ToString()) * Math.PI / 180;
                    LowY = Convert.ToDouble(pts.Rows[i - 1][1].ToString()) * Math.PI / 180;
                    MiddleX = Convert.ToDouble(pts.Rows[i][0].ToString()) * Math.PI / 180;
                    MiddleY = Convert.ToDouble(pts.Rows[i][1].ToString()) * Math.PI / 180;
                    HighX = Convert.ToDouble(pts.Rows[i + 1][0].ToString()) * Math.PI / 180;
                    HighY = Convert.ToDouble(pts.Rows[i + 1][1].ToString()) * Math.PI / 180;
                }
                AM = Math.Cos(MiddleY) * Math.Cos(MiddleX);
                BM = Math.Cos(MiddleY) * Math.Sin(MiddleX);
                CM = Math.Sin(MiddleY);
                AL = Math.Cos(LowY) * Math.Cos(LowX);
                BL = Math.Cos(LowY) * Math.Sin(LowX);
                CL = Math.Sin(LowY);
                AH = Math.Cos(HighY) * Math.Cos(HighX);
                BH = Math.Cos(HighY) * Math.Sin(HighX);
                CH = Math.Sin(HighY);
                CoefficientL = (AM * AM + BM * BM + CM * CM) / (AM * AL + BM * BL + CM * CL);
                CoefficientH = (AM * AM + BM * BM + CM * CM) / (AM * AH + BM * BH + CM * CH);
                ALtangent = CoefficientL * AL - AM;
                BLtangent = CoefficientL * BL - BM;
                CLtangent = CoefficientL * CL - CM;
                AHtangent = CoefficientH * AH - AM;
                BHtangent = CoefficientH * BH - BM;
                CHtangent = CoefficientH * CH - CM;
                AngleCos = (AHtangent * ALtangent + BHtangent * BLtangent + CHtangent * CLtangent) / (Math.Sqrt(AHtangent * AHtangent + BHtangent * BHtangent + CHtangent * CHtangent) * Math.Sqrt(ALtangent * ALtangent + BLtangent * BLtangent + CLtangent * CLtangent));
                AngleCos = Math.Acos(AngleCos);
                ANormalLine = BHtangent * CLtangent - CHtangent * BLtangent;
                BNormalLine = 0 - (AHtangent * CLtangent - CHtangent * ALtangent);
                CNormalLine = AHtangent * BLtangent - BHtangent * ALtangent;
                if (AM != 0)
                    OrientationValue = ANormalLine / AM;
                else if (BM != 0)
                    OrientationValue = BNormalLine / BM;
                else
                    OrientationValue = CNormalLine / CM;
                if (OrientationValue > 0)
                {
                    Sum1 += AngleCos;
                    Count1++;
                }
                else
                {
                    Sum2 += AngleCos;
                    Count2++;
                }
            }
            double tempSum1, tempSum2;
            tempSum1 = Sum1 + (2 * Math.PI * Count2 - Sum2);
            tempSum2 = (2 * Math.PI * Count1 - Sum1) + Sum2;
            if (Sum1 > Sum2)
            {
                if ((tempSum1 - (Count - 2) * Math.PI) < 1)
                    Sum = tempSum1;
                else
                    Sum = tempSum2;
            }
            else
            {
                if ((tempSum2 - (Count - 2) * Math.PI) < 1)
                    Sum = tempSum2;
                else
                    Sum = tempSum1;
            }
            totalArea = (Sum - (Count - 2) * Math.PI) * Radius * Radius;

            return totalArea; //返回总面积
        }
    }
}
